disp('flf');

x = 1.2;
kc = 2;

format long;
disp(fn_f2(x, kc));
disp(fn_f4(x, kc));

n = rand();
kc = rand();
k = rand();

disp(fn_int_tail(k, kc, n));
disp(fn_f2(-2*k/n, 2*kc/n)*8/2/pi/n^2);
disp(fn_f2(-2*k/n, 2*kc/n)*8/2/pi/n^2 - fn_int_tail(k,kc,n))

H = 0.01;
T = 0.01;
mu = 1.0;

h0 = H/T;

F = @(x, T) -T*log(1+exp(-x/T));
F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
    ./(1 + fn_inveta(m, h0)));


Nc = 10;

vx = linspace(0, 10, 400);
vep = vx.^2 - mu - T*log(2*cosh(H/2/T));
vep_ori = vep;

% please notice that the mu is double used
mu = zeros(Nc, length(vx));

while 1
    
    vep0 = vep;
    mu0 = mu;
    
    %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
    vep = vep_ori + fn_anf(2, vx, F(vep,T), vx);
    
    for n = 1:Nc
        
        vu = mu(n,:);
        vep = vep - T*fn_anf(n, vx, F2(vu,n), vx);
        
    end
    
    %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    for n = 1:Nc
        
        vu = -T*fn_anf(n, vx, F(vep, T), vx);
        for m = 1:Nc
            
            vu = vu + fn_Tnm(n,m,fn_inveta(m, h0),vx,mu(m,:),vx);
            
        end
        
        mu(n,:) = vu;
        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    fprintf('diss = %e \n', fn_norm(vep - vep0));
    fprintf('diss_m = %e \n', fn_norm(mu - mu0));
    
%     figure(1); hold on;
%     plot(fn_ex1(vx), fn_ex2(vep), '-o');
%     for n = 1:Nc
%         figure(n+1);
%         plot(fn_ex1(vx), fn_ex2(mu(n,:)), '-bo');
%     end
    
    if fn_norm(vep - vep0)<1e-8
        break;
    end
    
end

vx_nt = zeros(1,length(vx));
vx_len = vx(2:end) - vx(1:end-1);
vx_nt(1:end-1) = 1/2*vx_len;
vx_nt(2:end) = vx_nt(2:end) + 1/2*vx_len;

pressure = -2*sum(vx_nt.*F(vep,T))/2/pi;

fprintf('pressure = %f \n', pressure);

figure(1); hold on;
plot(fn_ex1(vx), fn_ex2(vep), '-o');
for n = 1:Nc
    figure(2); hold on;
    tv = mu(n,:) + fn_yn_inf(n, h0);
    fprintf('n = %d, y_inf =%f \n', n, fn_yn_inf(n,h0)); 
    plot(fn_ex1(vx), fn_ex2(tv), '-bo');
end

v_con_eta = zeros(1, Nc);
close(figure(4));
for n = 1:Nc
    
    v_str = mu(n, :);
    tv = fn_anf(n, vx, F2(v_str, n), vx);
    
    figure(4); hold on;
    plot(fn_ex1(vx), fn_ex2(tv), '-o');
    
    tmp = fn_int_with_tail(vx, tv);
    fprintf('n = %d, contribution = %f \n', n, tmp);
    v_con_eta(n) = tmp;
    
end

figure(64); hold on;
plot(1:Nc, v_con_eta, '-o');

Pn = fn_Pn(H, T, mu);
Qn = fn_Qn(H, T, mu);

close(figure(77));
close(figure(88));
close(figure(99));

figure(77);
hold on;
figure(88);
hold on;
figure(99);
hold on;

for n = 1:Nc
    
    figure(77);
    plot(fn_ex1(vx), fn_ex2(mu(n,:)), '-o');
    figure(88);
    plot(fn_ex1(vx), fn_ex2(Pn(n,:)), '-o');
    figure(99);
    plot(fn_ex1(vx), fn_ex2(Qn(n,:)), '-o');
    
end

% mu = 1.0;
% H = 0.01;
% T = 0.01;
% Nc = 4;
% disp(fn_pressure(mu, H, T, Nc));
% 
% mu = 1.0;
% H = 0.01;
% T = 0.01;
% Nc = 6;
% disp(fn_pressure(mu, H, T, Nc));
% 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function res = fn_G(x, T)

if x/T>33
    
    res = 0;
    
else

    res = 1./(1 + exp(x/T));
    
end

end



function Pn = fn_Pn(H, T, un)

[Nc, Num] = size(un);
Pn = zeros(Nc, Num);

if H>0

    for n = 1:Nc
        
        vu = un(n,:);
        P = (1 - exp(-vu))/T.*((n+1)*coth(H/2/T*(n+1)) - coth(H/2/T)) ...
            .*2.*sinh(H/2/T)^2 ...
            ./(cosh(H/2/T*(2*n+2)) - exp(-vu) - (exp(-vu) - 1).*cosh(H/2/T*2));
        Pn(n,:) = P;
        
    end

end

end

function Qn = fn_Qn(H, T, un)

[Nc, Num] = size(un);
Qn = zeros(Nc, Num);
h0 = H/T;

for n = 1:Nc
    vu = un(n,:);
    eta = fn_eta(n, h0);  
    Q = 1./(eta*exp(vu) + 1);
    Qn(n,:) = Q;
end

end



function pressure_dh = fn_pressure_dh(mu, H, T, Nc)

h0 = H/T;

F = @(x, T) -T*log(1+exp(-x/T));
F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
    ./(1 + fn_inveta(m, h0)));

vx = linspace(0, 10, 400);
vep = vx.^2 - mu - T*log(2*cosh(H/2/T));
vep_ori = vep;
% please notice that the mu is double used
mu = zeros(Nc, length(vx));
while 1
    
    vep0 = vep;
    mu0 = mu;
    
    %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
    vep = vep_ori + fn_anf(2, vx, F(vep,T), vx);   
    for n = 1:Nc     
        vu = mu(n,:);
        vep = vep - T*fn_anf(n, vx, F2(vu,n), vx);     
    end
    
    %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
    for n = 1:Nc
        vu = -T*fn_anf(n, vx, F(vep, T), vx);
        for m = 1:Nc       
            vu = vu + fn_Tnm(n,m,fn_inveta(m, h0),vx,mu(m,:),vx);       
        end
        mu(n,:) = vu;
    end
    
    fprintf('diss = %e \n', fn_norm(vep - vep0));
    fprintf('diss_m = %e \n', fn_norm(mu - mu0));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    if fn_norm(vep - vep0)<1e-8 && fn_norm(mu - mu0)<1e-8
        break;
    end
    
end

vx_nt = zeros(1,length(vx));
vx_len = vx(2:end) - vx(1:end-1);
vx_nt(1:end-1) = 1/2*vx_len;
vx_nt(2:end) = vx_nt(2:end) + 1/2*vx_len;

Pn = fn_Pn(H, T, mu);
Qn = fn_Qn(H, T, mu);

Num = length(vx);
vep_dh_ori = 1/2*tanh(H/2/T)*ones(1, Num);
vep_dh = vep_dh_ori;
mu_dh = zeros(Nc, Num);

while 1
    
    vep_dh0 = vep_dh;
    mu_dh0 = mu;
    
    %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
    vep_dh = vep_dh_ori + fn_anf(2, vx, fn_G(vep,T).*vep_dh, vx);   
    for n = 1:Nc     
        vep_dh = vep_dh - T*fn_anf(n, vx,Pn(n,:)-Qn(n,:).*mu_dh(n,:), vx);     
    end
    
    %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
    for n = 1:Nc
        vu = -T*fn_anf(n, vx, fn_G(vep,T).*vep_dh, vx);
        for m = 1:Nc       
            vu = vu + fn_Tnm0(n,m,vx,Pn(m,:)-Qn(m,:).*mu_dh(m,:),vx);       
        end
        mu_dh(n,:) = vu;
    end
    
    fprintf('diss = %e \n', fn_norm(vep_dh - vep_dh0));
    fprintf('diss_m = %e \n', fn_norm(mu_dh - mu_dh0));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    if fn_norm(vep_dh - vep_dh0)<1e-8 && fn_norm(mu_dh - mu_dh0)<1e-8
        break;
    end
    
end

pressure_dh = 2*sum(v_nt.*vep_dh.fn_G(vep));

end



function pressure = fn_pressure(mu, H, T, Nc)

h0 = H/T;

F = @(x, T) -T*log(1+exp(-x/T));
F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
    ./(1 + fn_inveta(m, h0)));

vx = linspace(0, 10, 400);
vep = vx.^2 - mu - T*log(2*cosh(H/2/T));
vep_ori = vep;
% please notice that the mu is double used
mu = zeros(Nc, length(vx));
while 1
    
    vep0 = vep;
    mu0 = mu;
    
    %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
    vep = vep_ori + fn_anf(2, vx, F(vep,T), vx);   
    for n = 1:Nc     
        vu = mu(n,:);
        vep = vep - T*fn_anf(n, vx, F2(vu,n), vx);     
    end
    
    %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
    for n = 1:Nc
        vu = -T*fn_anf(n, vx, F(vep, T), vx);
        for m = 1:Nc       
            vu = vu + fn_Tnm(n,m,fn_inveta(m, h0),vx,mu(m,:),vx);       
        end
        mu(n,:) = vu;
    end
    
    fprintf('diss = %e \n', fn_norm(vep - vep0));
    fprintf('diss_m = %e \n', fn_norm(mu - mu0));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    if fn_norm(vep - vep0)<1e-8 && fn_norm(mu - mu0)<1e-8
        break;
    end
    
end

vx_nt = zeros(1,length(vx));
vx_len = vx(2:end) - vx(1:end-1);
vx_nt(1:end-1) = 1/2*vx_len;
vx_nt(2:end) = vx_nt(2:end) + 1/2*vx_len;

pressure = -2*sum(vx_nt.*F(vep,T))/2/pi;

end


function res = fn_int_with_tail(vx, vf)

vx_nt = zeros(1,length(vx));
vx_len = vx(2:end) - vx(1:end-1);
vx_nt(1:end-1) = 1/2*vx_len;
vx_nt(2:end) = vx_nt(2:end) + 1/2*vx_len;

res =sum(vx_nt.*vf);

kc = vx(end);
v_tail = vx(end-9:end);
v2 = 1./v_tail.^2;
v4 = 1./v_tail.^4;
vf_tail = vf(end-9:end);
[alpha, beta] = fn_fit(v2, v4, vf_tail);

res = res + alpha/kc + beta/3/kc^3;

res = res*2;

end



function res = fn_yn_inf(n, h0)

if h0>0
    
    res = n*h0 ...
        + log((1+exp(-(2*n+2)*h0)-exp(-n*h0)-exp(-(n+2)*h0)) ...
        ./(1+exp(-2*h0)-2*exp(-h0)));

else
   
    res = log(n.^2 + 2*n);
    
end

end


function vu_out = fn_Tnm(n, m, a, vx, vu, vx_out)

vu_out = 0;

vf = log((1 + exp(-vu)*a)./(1 + a));

if m == n
    
    for i = 2:2:(2*n-2)
        
        vf_out = fn_anf(i, vx, vf, vx_out);
        vu_out = vu_out + 2*vf_out;
        
    end
    
    vu_out = vu_out + fn_anf(2*n, vx, vf, vx_out);
    
else
    
    vu_out = vu_out + fn_anf(abs(m-n), vx, vf, vx_out);
    vu_out = vu_out + fn_anf(m+n, vx, vf, vx_out);
    
    for i = (abs(m-n)+2):2:(m+n-2)
        
        vu_out = vu_out + 2*fn_anf(i, vx, vf, vx_out);
        
    end
    
end

end



function vu_out = fn_Tnm0(n, m, vx, vf, vx_out)

vu_out = 0;

if m == n
    
    for i = 2:2:(2*n-2)
        
        vf_out = fn_anf(i, vx, vf, vx_out);
        vu_out = vu_out + 2*vf_out;
        
    end
    
    vu_out = vu_out + fn_anf(2*n, vx, vf, vx_out);
    
else
    
    vu_out = vu_out + fn_anf(abs(m-n), vx, vf, vx_out);
    vu_out = vu_out + fn_anf(m+n, vx, vf, vx_out);
    
    for i = (abs(m-n)+2):2:(m+n-2)
        
        vu_out = vu_out + 2*fn_anf(i, vx, vf, vx_out);
        
    end
    
end

end




function res = fn_inveta(m, h0)

if h0>0
    res = sinh(h0/2)^2/(sinh(h0/2*(m+1))^2 - sinh(h0/2)^2);
else
    res = 1/(m^2 + 2*m);
end
     
end


function res = fn_eta(m, h0)

if h0>0
    res = (sinh(h0/2*(m+1))^2 - sinh(h0/2)^2)/sinh(h0/2)^2;
else
    res = m^2 + 2*m;
end

end



function vf_out4 = fn_anf(n, vx, vf, vx_out)

an = @(x) 1/2/pi*n./((n/2).^2 + x.^2);

vk = vx_out;
vp = vx;

mk = transpose(vk)*ones(1,length(vp));
mp = ones(length(vk),1)*vp;
ma = an(mk - mp) + an(mk + mp);

vpd = vp(2:end) - vp(1:end-1);
vp_nt = zeros(1,length(vp));
vp_nt(1:end-1) = 1/2*vpd;
vp_nt(2:end) = vp_nt(2:end) + 1/2*vpd;

vf_out = ma*diag(vp_nt)*transpose(vf);
vf_out = transpose(vf_out);
vf_out0 = vf_out;

kc = vp(end);
v_tail = vp(end-9:end);
v2 = 1./v_tail.^2;
v4 = 1./v_tail.^4;
vf_tail = vf(end-9:end);
[alpha, beta] = fn_fit(v2, v4, vf_tail);

%  fprintf('alpha = %f, beta = %f \n', alpha, beta);

vf_out_addon2 = alpha ...
    *8/2/pi/n^2*(fn_f2(2*vk/n, 2*kc/n) + fn_f2(-2*vk/n, 2*kc/n));

vf_out_addon4 = beta ...
    *32/2/pi/n^4*(fn_f4(2*vk/n, 2*kc/n) + fn_f4(-2*vk/n, 2*kc/n));

%  vf_out2 = vf_out0 + vf_out_addon2;

vf_out4 = vf_out0 + vf_out_addon2 + vf_out_addon4;

end


function vx_ex1 = fn_ex1(vx)

vx_ex1 = [-vx(end:-1:2), vx];

end

function vx_ex2 = fn_ex2(vx)

vx_ex2 = [vx(end:-1:2), vx];

end


function res = fn_int_tail(k, kc, n)

res = 8/pi*n*log(4)*k./(4*k.^2 + n^2).^2 ...
    + (4*k.^2 - n^2)*2./(4*k.^2 + n^2).^2 ...
    + 2/pi*n./(4*k.^2 + n^2)./kc ...
    - 4/pi*(4*k.^2 - n^2)./(4*k.^2 + n^2).^2.*atan(2*(k + kc)./n) ...
    + 8/pi*n*k./(4*k.^2 + n^2).^2.*log(kc.^2./(n^2 + 4*(k + kc).^2));

end


function res = fn_f2(x, kc)

res = -pi/2*(1 - x.^2)./(1 + x.^2).^2 ...
    + 1/kc./(1 + x.^2) + (1 - x.^2)./(1 + x.^2).^2.*atan(kc - x) ...
    + x./(1 + x.^2).^2.*log((1 + (kc - x).^2)./kc^2);

end


function res = fn_f4(x, kc)

res = 3*pi*(1 - 6*x.^2 + x.^4) ...
    + 6*(1 - 6*x.^2 + x.^4).*atan(x - kc) ...
    + 12*x.*(x.^2 - 1).*log((1 + (x - kc).^2)./kc^2) ...
    + 2*(1 + x.^2).^3/kc^3 ...
    + 6*x.*(1 + x.^2).^2/kc^2 ...
    + 6*(1 + x.^2).*(3*x.^2 - 1)/kc;

res = res./6./(1 + x.^2).^4;

end



function [alpha, beta] = fn_fit(v2, v4, v)

m = [v2*transpose(v2), v2*transpose(v4); ...
    v4*transpose(v2), v4*transpose(v4)];

RHS = [v2*transpose(v); v4*transpose(v)];

tv = m\RHS;
alpha = tv(1);
beta = tv(2);

end


function res = fn_norm(vx)

res = mean(abs(reshape(vx, [1, numel(vx)])));

end
